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The seismological dynamics of magnetars is largely determined by a strong hydro-magnetic coupling between the sohd 
crust and the fluid core. In this paper we set up a "spectral" computational framework in which the magnetar' s motion 
is decomposed into a series of basis functions which are associated with the crust and core vibrational eigenmodes. A 
general-relativistic formahsm is presented for evaluation of the core Alfven modes in the magnetic-flux coordinates, 
as well for eigenmode computation of a strongly magnetized crust of finite thickness. By considering coupling of 
the crustal modes to the continuum of Alfven modes in the core, we construct a fully relativistic dynamical model 
of the magnetar which allows: i) Fast and long simulations without numerical dissipation, ii) Very fine samphng of 
the stellar structure. We find that flie presence of strong magnetic field in the crust results in localizing of some high- 
frequency crustal elasto-magnetic modes with the radial number r?, > 1 to the regions of the crust where the field is 
nearly horizontal. While the hydro-magnetic coupling of these localized modes to the Alfven continuum in the core is 
reduced, their energy is drained on a time-scale of <C 1 s. Therefore flie puzzle of QPOs with frequencies larger than 
600 Hz stiU stands. 



Subject headings: Neutron stars 



1. Introduction 



> 

o 

(N 

d 



Magnetar oscillations have been subject of extensive theoret- 
ical research since the discovery of quasi-periodic oscillations 
(QPOs) in the hght curves of giant flares from soft gamma re- 
peaters (SGR) (Israel et al. 2005; Strohmayer & Watts 2005; 
Watts & Strohmayer 2006; see also Barat et al. 1983). The 
observed oscillations are measured with high signal-to-noise 
ratios during time intervals of typically few minutes in the fre- 
quency range between 18 and 1800 Hz. It has been proposed by 
many authors that the physical origin of the QPOs are seismic 
vibrations of the star; an idea which opens the possibiUty to 
perform asteroseismological analysis of neutron stars, giving a 
unique observational window into the stellar interior. Initially it 
was hypothesized that the observed oscillations originate from 
torsional shear modes which are confined in the magnetar crust 
(e.g. Duncan 1998, Piro 2005; Watts & Strohmayer 2006; 
Samuelsson & Andersson 2007; Watts & Reddy 2007; Steiner 
& Watts 2009). If this hypothesis were true, then the observed 
QPOs would strongly constrain physical parameters in the neu- 
tron star crust. However, it was soon reaUzed that, due to flie 
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presence of ultra strong magnetic fields (B ^ 10^* — 10^^ G; 
Kouveliotou et al., 1999) which are frozen both in the crust and 
the core of the star, the crustal motion is strongly coupled to 
the fluid core on timescales <C 1 s Levin (2006, hereafter L06). 
Over the years several authors have studied the coupled crust- 
core problem (Glampedakis, Samuelsson & Andersson 2006; 
Levin 2007, hereafter L07; Gruzinov 2008; Lee 2008; van 
Hoven & Levin 2011, hereafter vHLll; Gabler et al. 2011a; 
Colaiuda & Kokkotas 2011; Gabler et al. 2011b). In particular 
L06 and L07 argued that for sufficiently simple magnetic field 
configurations (i.e. axisymmetric poloidal fields), the Alfven- 
type motions on different flux surfaces are decoupled so that the 
Alfven frequencies in the core feature a continuum. This result 
is weU known from previous magnetohydrodynamic (MHD) 
studies, and it applies to general axisymmetric poloidal-toroidal 
magnetic fields (Poedts et al. 1985). It allows one to describe 
the problem of magnetar dynamics in terms of discrete crustal 
modes that couple to a continuum of Alfven modes in the core. 
With this approach, L07 and vHLll demonstrated that the pres- 
ence of an Alfven continuum has some important imphcations 
for magnetar oscillations: (i) Global modes of the star with fre- 
quencies that are located inside the continuum undergo strong 
exponential damping (this phenomenon is often called reso- 
nant absorption in the context of MHD (Goedbloed & Poedts 
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2004)). (ii) After the initial period (< 1 s) of exponential decay, 
the system tends to settle in a steady state in which it oscil- 
lates at frequencies close to the edges of the continuum; these 
oscillations correspond to the so-called edge-modes, that were 
first seen numerically in L07 and Gruzinov 2008, and were 
explained analytically in vHll. The edge-modes were further 
observed in the simulations of Gabler et al. (2011a) Colaiuda 
& Kokkotas (2011) and Gabler et al. (201 lb). 

In the past half-decade, two distinct computational strategies 
have been applied to the problem of calculating magnetar oscil- 
lations. 

(1) Several groups employed general relativistic MHD grid 
codes to simulate the dynamics of magnetized neutron stars. 
Sotani et al. 2008; Colaiuda et al. 2009 and Cerda-Duran et 
al. 2009 were able to reproduce continuum Alfven modes in 
the purely fluid stars with axisymmetric poloidal magnetic field, 
which provided important benchmark tests for the ability of the 
codes to handle complex MHD oscillations. Building on this, 
Gabler et al. (201 la), Colaiuda & Kokkotas (201 1) and Gabler 
et al. (2011b) included a crust in their neutron star models, and 
were thus able to study the coupled dynamics of the crust and 
the core. (2) Our group (L07 and vHLll) and Lee (2008) de- 
composed the motion of a magnetar into a set of basis functions, 
and studied the dynamics of the coefficients of these series ex- 
pansion; we shall refer to this strategy as the "spectral method". 
This framework is able to handle both the dynamical simula- 
tions and the stationary eigenmode problem; the latter reduces 
to solving the eigenvalue problem for a large matrix. L07 and 
vHl 1 chose the basis functions so that the crustal motion is de- 
composed into the normal modes of the free crust, and the core 
motion is decomposed into the sum of core Alfven modes and a 
separate contribution of the core's "dc" displacements in reac- 
tion to the motion of the crust. We refer the reader to Sections 
3.2 of L07, 4.2 of vHLl 1, and 4.2 of this paper for technical de- 
tails. This choice of basis functions casts the dynamics of mag- 
netars as a problem of coupled harmonic oscillators, in which 
the discrete modes of the crust are coupled to the Alfven modes 
in the core. 

The computations of vHl 1 have been performed using New- 
tonian equations of motion and in the limit of a thin crust. In 
this paper we improve on vHLll in two ways: 1) We adapt a 
realistic crust of finite thickness, threaded with a strong mag- 
netic field. 2) We employ fully relativistic equations govern- 
ing the motion of axial perturbations in the crust and the core. 
Our spectral method has several practical and conceptual advan- 
tages: (i) it is numerically inexpensive, making long simulations 
of the magnetar dynamics implemented on an ordinary worksta- 
tion possible, (ii) It allows one to sample the stellar structure at 
high spatial resolution, (iii) It does not suffer from the prob- 



lem of numerical viscosity that occurs in some finite difference 
schemes (scaling with the grid size), and it is able to handle ar- 
bitrary axisymmetric poloidal fields, and not just those that are 
the solutions of the Grad-Shafranov Equation^ 

The paper's plan is as follows. In section 2 we derive rela- 
tivistic equations describing the magnetic forces acting on ax- 
ial perturbations inside a neutron star with an axi-symmetric 
poloidal magnetic field. We construct a coordinate system 
which has one of its axes parallel to the fieldlines. The equa- 
tions thus obtained will in later sections when we calculate 
elasto-magnetic modes of the crust, and when we calculate the 
Alfven continuum in the core. 

In section 3.1 we introduce a formalism which allows us to 
calculate general relativistic elasto-magnetic eigenmodes of the 
crust by expanding the elasto-magnetic equations of motion in 
a set of basisfunctions. This reduces the eigenmode problem of 
the crust to a matrix eigenvalue problem. In sections 3.2 and 3.3 
we work out the relativistic equations describing the magnetic 
and elastic restoring-force densities in the curved space-time of 
the neutron star crust. In section 3.4 we apply these equations to 
the formalism of section 3.1 in order to find free crustal eigen- 
modes and -frequencies. 

In section 4, we find the core continuum Alfven modes in full 
general relativity, and we calculate their coupling to the crustal 
modes of section 3. The magnetar model constructed in this 
way, qualitatively shows the same features of the vHLl 1 model, 

1. e. above the fundamental Alfven frequency of ^ 20 Hz, the 
frequency domain is covered by the core continuum which ef- 
fectively acts to damp crustal motion. For particular choices of 
the field configuration, the continuum may contain a number 
of gaps, generally well below 200 Hz. These gaps give rise 
to the characteristic 'edge-modes' of vHLll. Moreover, the 
crustal modes that reside inside gaps remain undamped. In the 
appendix we revisit the problem of crustal mode damping due 
to the presence of an Alfven continuum, by analytically calcu- 
lating damping rates according to Fermi's golden rule. 

2. Relativistic equations for magnetic forces 

Magnetic coordinates 
We shall consider strongly sub-equipartition B <C lO^^G mag- 
netic fields, so that the physical deformation of the star is very 
small and the space-time is spherically-symmetric with respect 



^The approach developed by Sotani et al (2008) and used in Colaiuda et al. 
(2009, 2011) casts the MHD equations in the core into a particularly simple 
form; see section 4.4 of Sotani et al. (2008). This transformation is possible if 
the poloidal field is the solution of the Grad-Shafranov (GS) equation. There 
is, however, no compelling reason why the GS equation should hold, since neu- 
tron stars feature very strong stable stratification due to the radial gradients in 
proton-to-neutron ratios (Goldreich & Reisenegger 1992, Mastrano et al., 2011) 
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to the star's center. The metric can be written in the standard 
Schwarzschild-type coordinates r, 6 and 0. It is natural, in 
analogy with the Newtonian treatments, to introduce the flux 
coordinate system in which one of the axes is parallel to the 
magnetic field lines (the precise meaning of this construction 
in relativity is described below). In the axisymmetric poloidal 
field geometry the magnetic field lines are located in planes of 
constant azimuthal angle (f>, which allows us to define the two 
'magnetic' coordinates x{i\ ^) ■i/;(r, 9), such that the (co- 
variant) vectors = d/dcj) and = d/dx are orthogonal to 
~ d/dip. In the flux coordinate system the metric is given 

by 



-gttdf + gxxdx^ + 94>i>dip'^ 



while the magnetic-field vector is given by 

B = B^e^. 

Here B is the 4-vector whose components are given by 



(1) 



(2) 



(3) 



and Vi, is the 4-velocity vector which for the stationary star is 
given by vt = guv* = y/-gtt, Vi = 0. 

Clearly, gu and ij^^ are identical to the corresponding 
Schwarzschild metric terms. 



9tt = 1 



2m(r) 



sin^ I 



(4) 



MHD approximation, and Gauss' law for magnetic fields, i.e. 
{v^B^ — v^B^).^ — 0. For a static equillibrium, i.e. vt = 
\J—gu and = (where the index i runs over the spatial in- 
dices). Gauss' law can be expressed in the more famiUar form 



B\, 



V~9 







(8) 



where g = det {gij)/gtt- This expression provides basis for a 
convenient map between magnetic fields of Newtonian and rel- 
ativistic stars. In the Newtonian case, the flux coordinates x 
and ip are functions of r and 9; we keep this functional form 
for the relativistic versions of x and ip. The expression in Eq 
(|8]l is valid both in the curved space-time and in the flat Eu- 
clidean space (with g^j replaced by the Euclidean metric terms) 
of the Newtonian star We can therefore use Eq ([8]l to convert 
the values of the Euclidean field. Be, the correct values of 
the magnetic field in curved space-time, Bs (the subscript E 
stands again for Euclidean, S for Schwarzschild): Eq. (IHll gives 



0. We thus obtain 



gs 



(9) 



which results in the relativistic poloidal magnetic field which 
is tangent to the flux surfaces = const and which satisfies 
the Gauss' law. (In the following we will drop the subscript 
S.) In this work, for concreteness, we take use the Newtonian 
configuration of the magnetic field generated by a current loop 
inside the neutron star and discussed in detail in vHLll. Other 
Newtonian configurations are readily mapped onto the relativis- 
tic configurations using the procedure that is specified above. 



Maxwell's equations 

The evolution of the magnetic field is described by Maxwell's 
equations. In curved space-time these read 



(5) 



In the ideal MHD limit, the elective field Ef^ = v^F^i, van- 
ishes so that the only contribution to the electromagnetic tensor 
comes from the magnetic field: 



(6) 



After some manipulation, the relations (|5]l and (|6]) yield the 
MHD equations for the magnetic field: 



{v^'B'' -v''B^')^^Q. 



(7) 



This equation entails both magnetic induction, which describes 
the flux freezing that characterizes magnetic fields in the ideal 



Euler equations 

The equations of motion are obtained by enforcing conservation 
of momentum, i.e. by projecting the conservation of energy- 
momentum 4-vector on the space normal to the 4-velocity 



where the projection tensor is given by 



(10) 



(11) 



T^'' is the stress-energy tensor for a magnetized fluid in the 
ideal MHD approximation, and can be expressed as 



T^"" = [p + P+ v'^v'' +[P 



B' 



B^'B" 
An 



(12) 



Here, p and P are the mass-density and pressure and B^ 
B^Bfj^ is the square of the magnetic field, where 
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the covariant component of the Lorentz in- 
variant magnetic field 4-vector {e^^\a is the four dimensional 
Levi-Civita symbol and F^'^ is the electromagnetic tensor). The 
equations of motion become 



and 



B dC^ 



(16) 



P 



/i^^ ( P 



An 

Here we have used the relation ti^w" 



These equations can be combined into a single one. After restor- 
ing a factor of c^, we find 



B'^B^ 
47r 



K I ^— I ^ (13) 



P 
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g^ivv'^v" = -1. Eq. 
([TSjl together with equation (|7]l provides a full description of 
(incompressible) motion of the magnetized fluid in a neutron 
star. 

Perturbation equations 

We are now ready to derive equations that describe the lin- 
earized motion of a small Lagrangian fluid displacement 
about the static background equillibrium of the star The per- 
turbed components of the velocity and the magnetic field 4- 
vectors, v'^^^^ and B^^^ are 



9tt 



B^ \ 
B d 



9 XX 



Qtt 
9xx 



d_ 
dx 



(17) 



where ^ — y/g^C^ is the physical displacement (in the 0- 
direction) in unit length. This equation describes Alfven waves, 
traveling along magnetic field lines in the curved space-time of 
a magnetar We checked that in the non-relativistic limit Eq. 
([TTjl reduces to the correct expression for Alfven waves in self- 
gravitating magnetostatic equilhbria (Poedts et al., 1985). 



''pert 
^pcrt 



^SBf" 



(14) 



where the first terms on the right hand side denote the unper- 
turbed equillibrium quantities, and the second terms on the 
right hand side denote the Eulerian perturbations associated 
with the displacement In our 'magnetic' coordinates the 
only non-zero component of the unperturbed magnetic field is 
B^ — B / yjQxx^ ^'^'l because the equillibrium star is static and 
non-rotating the only non-zero component of the 4-velocity is 
= l/V— .gt*. Restricting ourselves to axi-symmetric tor- 
sional oscillations of the star, we introduce a small incom- 
pressible axisymmetric displacements C"^. This implies that 
Upg^^ = iv^':^^ = &v^-t, and that the perturbations in pressure 
5P and mass-density bp vanish. Technically, a full descrip- 
tion of the linearized motion of a neutron star would involve 
perturbations of the metric g^^y, requiring one to augment the 
above equations of motion with the perturbed Einstein equa- 
tions. However, since we're considering incompressional axial 
oscillations only, the metric perturbations are dominated by the 
current dipole moment. One can show that this causes pertur- 
bations in the off -diagonal elements of the metric tensor which 
are of order Sv^, so that the metric perturbations can be safely 
ignored (the so-called Cowling approximation). Taking these 
considerations into account, we linearize Eq's ( [T3] l and (|7]l and 
after some work we obtain 



P 



IP 

An 



9V 



9tt 



B d 



9xx '^'^9h 



3. Modes of a magnetized crust in General Relativity 

In this section we will describe a formalism that allows us 
to calculate relativistic eigenmodes and -frequencies of a neu- 
tron star crust of finite thickness and realistic equation of state, 
threaded with an arbitrary magnetic field. By considering a 
crust of finite thickness, we will obtain high frequency radial 
harmonics that are not present in the crust model of vHLl 1 but 
which should be taken into account in view of the observed high 
frequency QPO's. In the past several authors carried out theo- 
retical analyses of torsional oscillations of neutron stars with 
a magnetized crust. Piro (2005), Glampedakis et al. (2006) 
and Steiner & Watts (2009) considered horizontal shear waves 
in a plane-parallel crust threaded by a vertical magnetic field, 
whereas Sotani et al. (2008), Gabler et al. (2011a), Colaiuda 
& Kokkotas (2011) and Gabler et al. (2011b), performed grid- 
based simulations of spherical, relativistic stars with dipole 
magnetic fields. Lee (2008) on the other hand, studied the 
Newtonian dynamics of spherical magnetic neutron stars, by 
decomposing the perturbed quantities into a set of basis func- 
tions, and following the dynamics of the expansion coefficients. 
Here we follow a strategy which is closely related to that of 
Lee (2008). In this section, we consider normal modes of the 
'free' magnetized neutron star crust, i.e. in the absence of 
external forces. The idea in this section is to decompose the 
perturbed quantities into a set of orthogonal basis functions. 
By substituting this expansion in the equation of motion, we 
obtain equations for the evolution of the expansion coefficients. 
The solution of the crustal eigenmode problem, are in this way 
{9<k4>^-9tt5B'f') (15:teduced to a matrix eigenvalue problem. The hydromagnetic 
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coupling of the crust normal modes obtained in this section, to 
the core Alfven modes, will be discussed in section 4. 

Formalism for finding crustal eigenmodes 

In a magnetized and elastic crust, the motion of a small torsional 
Lagrangian displacement away from equillibrium ^{x,t) (we 
use the notation from vHLll; ^ denote crustal displacements, 
^ denote displacements in the core), can be described in the 
general form 



^mag 



(0 



(18) 



where iei and Lmag are the accelerations due to the elastic 
and magnetic forces acting on the displacement field. Expres- 
sions for and Lamg are given and discussed in the next 
sub-section. Augmented with no-tangential-stress conditions 
STriji = 6Tr0 = on the inner- and outer boundaries, this equa- 
tion describes the free oscillations of a magnetized neutron star 
crust. Our procedure for solving Eq. ( fTSj l is as follows: 

First, we decompose the crustal displacement field £,{t,x) 
into a set of basis functions 'i>i{x). 



(19) 



The functions form an orthonormal basis for a Hilbert space 
with inner product 



('7 10 - / w{x)rJ-Cd^x 
Jv 



(20) 



where ff and ( are arbitrary functions defined in the volume V 
of the crust, and w{x) is a weight function. Orthonormality of 
'ifi{x) implies that (^1/^ | = Sij, where Sij is the Kronecker 



delta. The coefficients aj of the expansion of Eq. ( 19 1 are then 
simply a,{t) = (|(t, f ) | ^^{x)). 

The next step is decompose the acceleration field of Eq. ( [TS] ) 



into basis functions according to Eq. ( 19 1, and to calculate 

the matrix elements {d^(_/dP \ ^j). This yields equations of 
motion for aj(i): 

'dj = Mij Qi, (21) 

where the double dot denotes double differentiation with respect 
to time, and where 



Clearly, a crustal eigenmode with frequency uj„i (i.e. a,„.j cx 
giw„t ^Qj. ^jj j-)^ jg jjQ^ simply an eigenvector of the matrix AI 



with eigenvalue — 



(22) 



The index m is used to label the different solutions to the above 
equation. In practical calculations, one truncates the series of 
Eq. ( 19 1 at a finite index i = N, so that one obtains a total 



number of N eigensolutions. The eigenvalue problem of Eq. 



(21 1 with finite (TV x TV) matrix M can be solved by means of 
standard linear algebra methods. Given a set of suitable basis 
functions, the eigenvectors and eigenvalues (or crustal eigen- 
frequencies) converge to the correct solutions of Eq. ( fTS] ! for 
sufficiently large iV (see the discussion of section 3.5). 

Orthogonality relation for elasto-magnetic modes 

In the limit of cx), the elasto-magnetic eigenfunctions are 



£,m{x) 



arn,t'^t{x), 



(23) 



where we omitted the time-dependent part e"^'"*, on both sides. 

The eigenfunctions will form a new basis for a Hilbert space 
of crustal displacements. We can introduce an inner product 
(...|...)i„o in which this basis is orthogonal as follows: Consider 
a deformation ^{x,t) of the crust, decomposed into a sum of 
eigenfunctions 



£,{x,t) = ^b,n{t)£,rn{x), 



(24) 



where we incorporated the harmonic time dependence in the 
coefficients bm{t). Since are the eigenmodes of the crust, 
the kinetic energy of the displacement field K{£^) must be equal 
to the sum of kinetic energies of the individual modes K{bm£,m) 

In the static Schwarzschild space-time of the neutron star, the 
conjugate time-like momentum pt — ~E is a. constant of 
geodesic motion (see e.g. Misner, Thorne & Wheeler (1973), 
§25.2). In terms of the locally measured energy i?L = \J—9ttV^, 
the conserved "redshifted" energy is £' —pt = \J—gttE\j. 
Similarly, the kinetic energy K in terms of the locally measured 
kinetic energy Kj^ is 



K 



t) = V^tKL (I) 



V-9ttP 



dr 



dV = 



dt 



dV (26) 



-M/dt\d£,/dt), 
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where p = [p + P/c^ + B'^/inc'^) is the mass-density in a 
local Lorentz frame, and dV ~ ^grrQAAQBR dr d4> dO is the lo- 



cally measured space-like volume element. By substituting this 



expression for the kinetic energy into Eq. (25 i, one finds that 
the cross-terms, {d^m/dt \ d£,k/dt)mc = i^m^k{£.m \ &)me 
with m k, vanish. After normalizing the eigenfunctions f™, 

so that Kibmira) 

lation: 



l/2ci;^„6^, we obtain the orthogonality re- 



(e™ I e 



kj mc 



p 



im-ikdV = 5, 



mk ■ 



(27) 



The coefficients 5,„(t) are now simply obtained by taking the 
inner product between the displacement field ^{x,t) and the 
eigenfunctions ^m{x)'. 



hm{t) = {i{x,t)\U{x)), 



(28) 



In the next two sections we give expressions for L^ag and 
Loi, and we discuss our choice of basis functions '^i and the 
resulting boundary forces (due to the no-stress boundary con- 
ditions) at the end of section (3.2). In section (3.3) we set up 
a realistic model of the magnetar crust and we calculate the 
corresponding elasto-magnetic modes in section (3.4), where 
we apply the formalism described above. In the remainder of 
this paper, we will focus solely on axi-symmetric azimuthal 
displacement fields, i.e. ^ = ^ (where is the unit vector 
in the azimuthal dkection and f is the displacement amplitude) 
and dl/dcj) = 0. 



where Ti„ag(r) is the magnetic stress at r, and e is an infinitesi- 
mal number. Adding the boundary terms, we obtain 



L 



mag(0 — 



B 



d 



9xx 47rc2p^5^ dx 



9tt „ d 
9xx 9x 




(30) 



where the 5's are Dirac delta functions. The magnetic stress 
Tmag is derived by linearizing Eq. ([T2ji and retaining first order 
terms. One obtains 



^/9tt9<p<P B 
cos a 



d 



9 XX 



47r dx V ^500 / 



(31) 



3.2. Relativistic equations for elastic forces 



In the following we use relativistic equations describing the 
elastic force density acting on axial perturbations in the crust 
as derived by Schumaker & Thorne (1983) (see also Karlovini 
& Samuelsson 2007), and presented in a convenient form by 
Samuelsson & Andersson (2007, SA) (for more details on the 
derivation of the following equations we refer the reader to these 
two papers). As shown in SA, the equation of motion for axial 
perturbations in a purely elastic crust, i.e. d^£,/dt^ — ici(0' 
can be solved by expanding the displacement field ^(r, 9, cj)) into 

vector spherical harmonics ^-ajmi^, (f) oc rx VY[^ (where Y™ 
is a spherical harmonic of degree I and order m) and corre- 
sponding radial- and time-dependent parts ^R(r) and frit) of 
the displacement field. Rewriting Eq. (2) of SA now gives 



3.1. Magnetic force density in the free crust 

While the equations of section 2 hold at arbitrary locations 
in the star, we will now consider magnetic forces acting on axi- 
symmetric, azimuthal perturbations £,{r,9) = £,{r,6)e^ in the 
'free' crust, i.e. a crust with no external stresses acting on it. 
This implies that to Eq. (17 1 we have to add boundary force 



terms arising from this no-external-stress condition. The tan- 
gential forces per unit area on both boundaries are given by 



■^mag 



(^n + e) - Tmag(nn - e) 



T;,ag(nn + e) (29) 



^mag(^out ^) ^mag(^o 



— ^mag ( 



mag V ' out 




gtt d 

9rr dr 

-mt 5 




(32) 



where the metric terms gtt and the standard Schwarz- 

schild metric terms, and /i(r) is the (isotropic) shear modulus. 

The expansion of ^ into vector spherical harmonics, leads to a 
particularly simple stress-free boundary condition for the radial 
function ^r : 



dr \ r 







which is valid on the inner- and outer boundaries, r 
r = ri. 



(33) 
ro and 
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We are now ready to select our basis functions in order 
to solve Eq. ( 18 1. It is convenient to seperate 'i'i into angular 
and radial parts, i.e. '^i = ^h.; "^r,!- Although our particu- 
lar choice of basis is technically arbitrary, in view of the above 
discussion a natural choice for the angular part ^h./ are vector 
spherical harmonics of order rn = and I = 2,4, 6... etc. (we 
consider axi-symmetric perturbations which are anti-symmetric 
with respect to the equator). 



47r 



' e, (34) 



which are orthonormal with respect to the following inner prod- 
uct: 



H,/ 



*H,/-*H,/' sin 61^61 = Su, 



(35) 



One tempting choice for the radial function is to use the radial 
eigenmodes of Eq. (33 i, ^r_„, (where n is the number of radial 
nodes) as basis functions, i.e. vPr^h = Cr.h- It turns out how- 
ever, that the expansion of the elasto-magnetic displacement 



field [see Eq. (19i] into elastic eigenfunctions is very ineffi- 



cient. We found that better convergence is realized with 



ri - ro 



ri - ro 



ri - ro 



for n=0 



(36) 



which obey Eq. ( 33 1, so that no extra boundary terms in Lc\ are 
needed to preserve the stress-free condition. The basis functions 
of Eq. ( [36| are orthonormal with respect to the following inner- 
product: 



(«'R,« I ^'R,«') = / *R,n *R,«'^dr Snn' 
J rn ^ 



(37) 



Combining Eq's (34i and (37i gives us a series of basis func- 
tions that we use in the next section to calculate elasto-magnetic 
modes of the crust 



(38) 



which are orthonormal 

/? , ? ■ Z"'' /""sine 



d9dr = SwSnn' (39) 



Note that the weight function w of Eq. ( 20 1 takes the form 

w{r, 9) — sinO/r^. 



3.3. The neutron star model 

We assume that our star is non-rotating and neglect defor- 
mations due to magnetic pressure, which are expected to be 
small. Therefore, we adopt a spherically symmetric background 
stellar model that is a solution of the Tolman-Oppenheimer- 
Volkoff equation (TOV equation). We calculate the hydro- 
static equillibrium using a SLy equation of state (Douchin & 
Haensel, 2001; Haensel & Potekhin, 2004; Haensel, Potekhin 
& Yakovlev, 2007) (see http://www.ioffe.ru/astro/NSG/NSEOS/ 
for a tabulated version). The model that we use throughout 
this paper has a mass of = 1.4 Mq, a radius i?* = 
1.16- 10^ cm, a crust thickness AR — 7.9- lO"* cm, a cen- 
tral density pc — 9.83- 10^"* g cm^"^ and cental pressure Pc = 
1.36- 10'^^ dyn cm~^. The crustal shear modulus /i is given by 
(Strohmayer et al., 1991) 



0.1194 



1 + 0.595(173/r)2 



(40) 



where n is the ion density, a = (3/47rn)^/^ is the average spa- 
cing between ions and F = [ZeY / aksT is the Coulomb cou- 
pling parameter. We evaluate /i in the limit F — > oo. 

To the spherical star we add a poloidal magnetic field, which 
we generate as follows: We start with an Euclidean (flat) space 
into which we place a circular current loop of radius rd = 
0.55 and current / and calculate the magnetic field gener- 
ated by the loop (see e.g. Jackson, 1998). Then we map this 
field onto the curved space-time of the neutron star, as discussed 
in section 2. The field is singular near the current loop, however 
all the field lines which connect to the crust (and thus are phys- 
ically related to observable oscillations) carry finite field val- 
ues. This particular field configuration is chosen as an example; 
there is an infinite number of ways to generate poloidal field 
configurations. In figure [T] we plot resulting shear- and Alfven 
velocities in the crust as a function of radial coordinate r. 

3.4. Results 

We now use the formalism and equations of the previous sec- 
tions to calculate elasto-magnetic modes of the magnetar crust. 
We construct a basis from Nn radial functions ^'r,„ (r) (see Eq. 
(j36|) with index n = 0, 1, Nn — 1, and Ni angular functions 
^'hj(6') (see Eq. Q) with even index ^ = 2, 4, 2Ni. These 



functions provide a set of Nn x Ni linearly independent basis- 
functions '^in. Using this basis, we solve the matrix equation 
(22 1, and reconstruct the normal modes according to Eq. ([T9]l. 



Radial and horizontal cross-sections of a selection of eigen- 
modes are plotted in figures |2] and |3] and table |2] contains a list 
of frequencies. These results are based on a stellar model with 
a poloidal field strength of lO^'^ G at the magnetic pole. For 
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2e+08 



1 .5e+08 




5e+07 



Fig. 1. — Shear velocity Cg — \f[ijp (solid line) versus Alfven ve- 
locity CA = v^BV^Trp for a poloidal field strength of 10^^ G (dotted 
line). The dashed lines are the radial components of the Alfven veloc- 
ity, CA,rad = CA COS ft, evaluated at (from left to right) Q = 69°, 79° 
and 89°. Closer to the poles (smaller &), the field becomes nearly ra- 
dial and CA.rad ~ CA. The CA-curve shown in this plot is evaluated at 
the pole {d — 0°), but varies negligibly as a function of 8. 



the calculation we used Nn = 35 radial basis functions and 
Ni = 35 angular basis functions. We labeled the modes with 
integer indices rti — 0, 1, 2... and li — 2, 4, 6, where ni is 
defined as the number of nodes along the r-axis and /i + 1 is the 
number of nodes along the 0-axis (including the poles). Note 
that the index li, in contrast to I, does not signify a spherical 
harmonic degree since the angular dependence of the elasto- 
magnetic modes differs from pure spherical harmonics. How- 
ever, there is a connection between the two indices: the elasto- 
magnetic mode of degree li and order rii, can be interpreted 
as the magnetically perturbed elastic mode of the same order 
and (spherical harmonic-) degree. More precisely, if one grad- 
ually increases the magnetic field strength, the n, I elastic mode 
transforms into the elasto-magnetic mode of the same indices, 
ni = n and li = I (see fig. |6]l. It is interesting to note (see fig's 
[6]and|3]l that as the field strength increases, modes with ?ii > 
become more and more confined to a narrow region near the 
equator (a similar effect was recently observed in the grid-based 
simulations of Gabler et al. 201 lb). In the equatorial regions, 
the horizontal field creates a magnetic tension-free cavity for 
modes with radial nodes, which are reflected back towards the 
equator at higher lattitudes where the field becomes more ra- 
diaj^ The ni = modes however, having no radial nodes. 



are virtually insensitive to the magnetic field and are therefore 
not confined to low lattitudes. The field strength-dependence of 
the eigenfrequencies is illustrated in figure |5] As we increase 
the field strength, we find that the increase in frequency 6uj for 
ni — modes scales weakly with B, i.e. Suj oc B^. For modes 
with Til > 0, (5w cx -B^ if B < 5- lO^^ q ^nd Suj (x B if 
B>5- 10^3 G. 



Table 1 : Normal mode frequencies 
mode indices elastic modes 
(B = G) 



elasto-magnetic modes 
(B = IQis G) 



ni 


= 0, 


h 


= 2 


27.42 Hz 


27.61 Hz 


ni 


= 0, 


h 


= 4 


58.16 Hz 


59.14 Hz 


ni 


= 0, 


h 


= 6 


86.69 Hz 


88.13 Hz 


ni 


= 0, 


h 


^ 8 


114.7 Hz 


116.5 Hz 


ni 


= 1, 


h 


= 2 


895.9 Hz 


954.1 Hz 


ni 


= 1, 


h 


= 4 


897.4 Hz 


985.7 Hz 


ni 


= 1, 


h 


= 6 


899.7 Hz 


1001.4 Hz 


ni 


= 1, 


h 


= 8 


902.8 Hz 


1003.4 Hz 


Til 


= 2, 


h 


= 2 


1474.6 Hz 


1607.1 Hz 


ni 


= 2, 


h 


= 4 


1475.7 Hz 


1664.4 Hz 


ni 


= 2, 


h 


= 6 


1477.5 Hz 


1708.1 Hz 


ni 


^ 2, 


h 


^ 8 


1479.9 Hz 


1740.4 Hz 



Table 2; The eigenfrequencies of the non-magnetic crust (second col- 
umn) versus the eigenfrequencies of the magnetized crust (third col- 
umn), with a magnetic field of 10^^ G at the polar surface. The elasto- 
magnetic frequencies were calculated using a basis of 35 x 35 basis- 
functions 'i'ln- 

As a test, we compared the eigenfrequencies and eigenmodes 
for zero field, i? = 0, to those obtained by a direct integra- 
tion of the elastic equation of motion of Eq. (33 i ^ We find 
that both frequencies and wavefunctions obtained by the series 
expansion-method converge rapidl}|^to real values, obtained by 
integration of Eq. ([33]l. E.g. for Nn — 10, rti = elastic 



similar effect is well-known from the study of waveguides: as tlie waveguide 
gets narrower (i.e. as its transverse frequency increases), the propagating wave 



may become evanescent in the longitudinal direction and be reflected 
•^Tlie latter works as follows: One starts by assuming harmonic time depence for 
the displacement so that icl(f) = — '^J'^f- Dropping the angular- and time- 
dependent parts of ^ on both sides of the equation, one is left with an equation 
for ^R^, which is integrated from the bottom of the crust, with corresponding 
boundary condition, to the surface. This is repeated for different uj until the sur- 
face boundary condition is satisfied; one has found an eigenmode. By repeating 
this procedure with gradually increasing ui, one obtains a series of eigenmodes 
and -frequencies. 

*Note that in the purely elastic case, I is a good quantum number and the angular 
basis functions ^h,;(^) are already solutions to the elastic eigenmode equation. 
Therefore, for a given Zi = I only the series with the radial basis-functions 
needs to be considered. 
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Fig. 2. — Radial profiles of h —2 elasto-magnetic modes, evaluated 
at9 — 81°. The vertical scale of individual curves is adapted for visual 
convenience. 
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Fig. 3. — Examples of elasto-magnetic eigenmodes for Bp = 10^^ G 
(where Bp is field strength at the magnetic pole), as a function of the 
polar angle 9, evaluated at the crust-core interface. The ni — modes 
are nearly unaffected by the magnetic field and are spread out over the 
crust, whereas the ni > modes are affected strongly by the magnetic 
field, and are confined to regions near the equator, where the field is 
horizontal. 
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ii=2,n,=1 
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11 11.1 11.2 11.3 11.4 11.5 11.6 
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Fig. 4. — Elastic crustal modes obtained through integration of the 
elastic equation of motion (thick dashed curves), and the same modes 
obtained by the series-expansion method (overplotted by the thin solid 
curve), using A''„ radial basis functions. 



frequencies have a typical error of 0.02%, while frequencies for 
modes rii < 4 are well within 1% accuracy. In figure |4] we 
plot elastic eigenfunctions, obtained by both methods. The so- 
lutions from the series-expansion method with Nn = 10 radial 
basis functions are nearly indistinguishable from the solutions 
obtained by direct integration. 



For the full elasto-magnetic equation of motion, Eq. (18 1 
with a magnetic field strength of 10^^ G at the pole, we tested 
the convergence of resulting eigenfrequencies by increasing the 
number of basis functions Nn and Ni (see figure [7|. We find 
that, compared to the non-magnetic case, a significant number 
Nn of radial functions and Ni angular functions is required to 
get acceptable convergence to stable results. The large num- 
ber of required radial basis functions can be understood from 



the fact that the magnetic acceleration L b (Eq. ( 3 1 1) contains 
delta-functions, arising from the boundary terms. Obviously, 
one needs many radial basis functions to obtain an acceptable 
sampling of these boundary terms. The number of computa- 
tional operations however, is a steep function of the number of 
basis functions (approximately cx {Ni x Nn)^), so that computa- 
tions with large Ni and N^ can become unpractical on ordinary 
workstations. Although this limits the number of basisfunctions 
in our calculations, we find that for Ni , Nn ~ 35, the scatter in 
frequencies is typically < 1% for most modes (figure [7|i, and 
the eigenfunctions reproduce the orthogonality relation of 



Eq. ( 27 1 with good precision. 
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Fig. 5. — Frequencies as a function of B. For ni > 0, the frequencies 
of (low) ii-modes nearly coincide and are therefore collectively indi- 
cated with their ni-value, i.e. ni = 1, ni = 2, etc. Note that high 
field strengths, the ni > frequencies collectively behave as a; oc B. 
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Fig. 6. — Angular geometry of the Zi = 2, ni = 1 crustal mode (at the 
crust-core interface), as a function of the magnetic field strength. For 
zero magnetic field, the curve is identical to the / — 2 vector spherical 
harmonic >I'h,;(^). As the field strength increases, the crustal motion 
becomes gradually more confined towards the equator 
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Fig. 7. — Demonstration of convergence for elasto-magnetic fre- 
quencies for low-order, low-degree modes as a function of Nn and 
Ni, where we took Nn = A^;. The actual number of basisfunctions, 
A'^ = Nn X Ni, is the square of the value along the x-axis. 

4. Core continuum and crust-core coupling 
4.1. The continuum 

The equation of motion is in this case simply the Alfven 
wave equation: 



= imag K(V',X)] 



(41) 



where t denotes the Schwarzschild time-coordinate. The oper- 
ator £mag is given in Eq. ([TTji, which we repeat here for conve- 
nience 




(42) 



Here gu, 51^^ and are the metric terms corresponding to the 
system of coordinates defined in section 2. 

For determining the spectrum of the core continuum, the ap- 
propriate boundary conditions are ^(x = Xc) = 0, where XcW 
marks the location of the crust-core interface. The full signifi- 
cance of this boundary condition will become apparent later in 
this section when we develop the analysis for the crust-core in- 
teraction; see also section 4.2 in vHLll. With this boundary 



condition. Equation (41 1 constitutes a Sturm-Liouville prob- 
lem on each separate flux surface t/j. Using the stellar struc- 
ture model and magnetic field configuration described in sec- 
tion 3.3, we can calculate the eigenfunctions and eigenfrequen- 
cies for each flux surface ip. The reflection symmetry of the 
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stellar model and the magnetic field with respect to the equa- 
torial plane assures that the eigenfunctions of equation ( [4T] i are 
either symmetric or anti-symmetric with respect to the equa- 
torial plane. We can therefore determine the eigenfunctions by 



integrating equation (41 1 along the magnetic field lines from the 
equatorial plane x = to the crust-core interface x — Xc (^)- 
Let us consider the odd modes here for which ^ (0) — 0, and 
solve equation (41 1 with the boundary condition ^ [xc) = at 



the crust-core interface; for even modes, the boundary condi- 
tion is (0) /dx = 0. We find the eigenfunctions by means of 
a shooting method; using fourth order Runge-Kutta integration 
we integrate from x = to x = Xc- The correct eigenvalues 
cr„ and eigenfunctions ^„ (x) are found by changing the value 
of a until the boundary condition at ^„ is satisfied. In this way 
we gradually increase the value of a until the desired number of 
harmonics is obtained. In figure |8] we show a typical resulting 
core-continuum. The continuum is piece-wise, and covers the 
domains a = [41.8,67.5] Hz and a [91.4, oo) Hz. Gaps, 
such as the one between 67.5 Hz and 91.4 Hz in fig. [8] are a 
characteristic feature for the type of poloidal field that we em- 
ploy in this paper, and typically occur at low frequencies (i.e. 
(7 < 150 Hz). As we discuss in section 4.3, they may give rise 
to strong low frequency QPOs; see also vHLl 1 and Colaiuda & 
Kokkotas 2011. 
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Fig. 8. — The curves show the Alfven frequencies (t„ as a function of 
the angle the polar angle at which the flux surface -0 intersects 

the crust. Since we are only considering odd crustal modes, the only 
Alfven modes that couple to the motion of the star are the ones with 
an odd harmonic number n. This particular continuum was calculated 
using a poloidal field with a surface value of B = 10^^ G at the poles. 

According to Sturm-Liouville theory the normalized eigen- 




functions ^„ of equation (41 1 form an orthonomial basis with 



respect to the following inner product: 

(Cm, in) = f r{x) Cm ix) Cn (x) dx = S,n,n (43) 

Jo 

Where 6m,n is the Kronecker delta. Noting that the operator 
imag(C) is in Sturm-Liouville form, one reads off the weight- 
function r(x): 



9xx 47rp 



(44) 



9tt 

We have checked that the solutions Cn(x) satisfy the ortho- 
gonality relations. 

4.2. Equations of motion for the coupled crust and core 

We are now ready to compute the coupled crust-core motion. 
In contrast to L07 and vHLl 1, where the crust was assumed to 
be an infinitely thin spherical elastic shell, we shall here adopt 
a crust of finite thickness with realistic structure. We label the 
lattitudinal location by the flux surface intersecting the crust- 
core interface, and consider the crustal axisymmetric displace- 
ments icjiii^, r), where r is the radial Schwarzschild-coordinate. 
In the MHD approximation, the magnetic stresses enforce a no- 
slip boundary condition at the crust-core interface (at r = rg 
in the Schwarzschild coordinates of the crust, or Xc in the flux- 
coordinates of the core), such that ^ (V': Xc) — £, (^("0): ^o) in- 
stead of C (V', Xc) = 0. It is useful to make the following sub- 
stitution 

C i^, x) = C (0, X) - HOW. ro) w x) (45) 

where we choose the function w{ip,x) so that (1) it corre- 
sponds to the static displacement in the core and hence satisfies 
F {w x)) = 0, and (2) w Xc) = 1- From the definition 
of the operator F it follows that for the odd modes 



' (V', x) 



9tt g^^B{'ip,x') 



(46) 



Here the constant K (ip) is chosen such that w (ip, Xc) = 1- 



The new quantity C from Eq. (45 i now satisfies the boundary 
condition {ip, Xc) = and can be expanded into the Alfven 
normal modes Cn which satisfy the same boundary conditions. 



We now proceed by substituting equation (45 i into equa- 
tion ( |4T] l thus obtaining a simple equation of motion for ( 



5\(V',x) 



Jmag 



(C(V',x)) = ~w{ip,x) 



We expand ( and w into a series of Cn's: 

C {fP, X,t) = ^ fln {-Ip, t) in {lp, X) 



w{ip,x) 



Cn (lp) in {lp, X) ■ 



(47) 

(48) 
(49) 
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Using these expansions, equation (47 1 reduces to the following 
equations of motion for the eigenmode amplitudes a„ 



d'^an (-0) 



(■0) a„ (0) 



(0) 



a<2 



(50) 



These equations show how the core Alfven modes are driven by 
the motion of the crust. To close the system, we must address 
the motion of the crust driven by the hydromagnetic pull from 
the core: 




cos a 



dx 



S{r ~ ro) 



(51) 



The expression between the square brackets is the hydro- 
magnetic stress from stellar core acting on the crust, a is the 
angle between the magnetic field line and the radial coordinate 
of the star and icrust (C) — -^mag (C) + -^ci (C) is the acceler- 
ation of the crustal displacement due to magnetic- and elastic 
stress (see section 3). We can rewrite this in terms of the coef- 



ficients, using Eq. (45 i, the definition of w, and the expansions 



and orthogonality relations of Eq's (27 1 and (29 1, as 



\/9rrgtt 



9xx 



2c2 



COS a 



E 



dx 



(52) 



9xx 



K 



9tt B. 



i'q sin 6d9 



where the coefficients bj (t) are crustal mode amplitudes defined 



in Eq's (24i and (29 1. Up to this point the derived equations of 



motion for the crust and the fluid core are exact. We are now 
ready to discretize the continuum by converting the integral of 
equation (jSTJi into a sum over N points 9i . In order to avoid the 
effect of phase coherence (see section 3) which caused drifts in 
the results of L07, we sample the continuum randomly over the 
6'-interval [0, 7r/2]. In the following, functional dependence of 
flie coordinate ip or 6 is substituted by the discrete index i 
which denotes the i-th flux surface. 



E 



+ = 

y/grr.i9tt,i Bf 



9xx, 



2c2 



COS a," 



E' 



' dx 



(53) 



/ 9xX,i 

9tta Bi^Jg^^ 



m / 



rn sin OjAOi 



3 



(54) 



These are the equations that fully describe dynamics of our 
magnetar model. As with the toy model from section 2 we in- 
tegrate them using a second order leap-frog scheme which con- 
serves the total energy to high precision. As a test we keep track 
of the total energy of the system during the simulations. Further 
we have checked our results by integrating equations (53) and 



( 54 1 with the fourth-order Runge-Kutta scheme for several runs 
and found good agreement with the leap-frog integration. 

4.3. Results 



Based on the results of vHLll, we expect the following dy- 
namical characteristics to occur; 1) Crustal modes with frequen- 
cies that are inside the continuum should undergo resonant ab- 
sorption, i.e. if such a mode couples efficiently to continuum 
Alfven modes of the core with similar frequencies, its motion 
will be damped on rather short time-scales. In the appendix 
we analytically investigate the efficiency of this coupling and 
the resulting damping time scales. 2) Late-time behavior of the 
system will show oscillations near the edges of the continuum; 
the edge modes. 3) Gaps, as present in the continuum of fig. 
[Hjwill give rise two types of QPOs. First, crustal modes which 
are inside these gaps will remain undamped, although slightly 
shifted in frequency due to the interaction with the continuum. 
Second, edge modes near the edges of the gaps may occur. 
All of these characteristics were observed in simulations of 
vHLll, and we expect them to occur in this work. 

We consider 16 crustal modes, i.e. (n, /) = (0,2), (0,4), 
(0,6), (0,8), (0,10), (0,12), (0,14), (0,16), (0,18), (0,20), 
(1, 2), (1, 4), (1, 6), (1, 8), (1, 10) and (1, 12). We couple these 
crustal modes to 9000 continuum oscillators, i.e. 300 different 
flux surfaces, each with 30 Alfven overtones. We start the simu- 
lation by initializing the crustal mode amplitudes bj = 1 for all 
crustal modes, while keeping the continuum oscillators relaxed 
{ttni = 0). We evolve the system for 52s in time. 
In figure |9] we show the power spectrum which was calculated 
using the data of the last 26s of the simulation. 

5. Discussion 

In this paper we have laid out the spectral formalism for com- 
putation of general-relativistic torsional magnetar oscillations. 
This method is efficient; a typical simulation of 50 seconds of 
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Fig. 9. — Power spectrum of the crustal motion. 
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Fig. 10. — The same power spectrum, close up. 
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Fig. 11. — Displacement of the h = 2, ni = 1 mode. The theoreti- 
cally calculated damping time is Td = 5.8- 10"'^ s. Note the transient 
increase in the mode amplitude. This is due to the initial Alfven wave 
train, which is reflected at the equator. 

the magnetar dynamics (i.e., up to tens of thousands of the os- 
cillatory periods) takes only a few hours or an ordinary worksta- 
tion. By contrast, published calculations by groups using alter- 
native finite-different schemes (Gabler et al. 201 la, Colaiuda & 
Kokkotas 201 1, Gabler et al. 201 lb) use massive computational 
resources but still can track at most several tens of the oscilla- 
tions periods. The second-order symplectic leap-frog scheme 
ensures that the energy of the system is conserved with very 
high accuracy. Our simulations allow us to investigate which 
of the oscillatory behavior is long-lived enough 100s) to be 
relevant to the observations of QPOS in the tails of giant SGR 
flares (Israel et al. 2005, Strohmayer & Watts 2006). 

One of the puzzling features of the observations are sev- 
eral high-frequency QPOs above 600 Hz (Watts & Strohmayer 
2006). The thin-crust models of vHll had strongly suggested 
that crustal modes of such high frequency should be subject to 
the strong resonant absorption in the core, even if the core's 
Alfven modes do not form a mathematical continuurr0 In ac- 
cordance with recent results of Gabler et al. (201 lb), we found 
that some crustal modes are confined to the regions in the crust 
where the magnetic field is nearly horizontal. Because of this, 
the coupling to the Alfven modes in the core is reduced relative 
to the coupling strength estimated in vHLl 1, however, the cou- 
pling is still large enough for the mode energy to be drained on 
a time-scale small compared to the observed QPOs (r^; ^ 100 



^This is because the frequencies of even discrete Alfven modes form a grid, 
whose characteristic spacing is much less than 600Hz. At such high frequen- 
cies, the grid acts dynamically as a continuum. See vHl 1 for a more detailed 
discussion 
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s). Thus the problem of the high frequency QPOs (> 600 Hz) Strohmayer, T. E. & Watts, A. L., 2005, ApJ, 632, Ll 1 1 

still stands. Watts, A. L. & Reddy, S. 2007, MNRAS Letters, 379, 63 

Watts, A. L. & Strohmayer T. E., 2006, ApJ, 637, L117 
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A. Damped modes 

Now we explore the phenomenon of resonant absorption which occurs in a system where a harmonic oscillator is coupled to a 
continuum of oscillators. Our aim is to find an analytic estimate for the rate at which the energy of such an oscillator is transferred 
to the continuum. The objective of this section and the method that we follow, are analogous to a derivation of the quantum 
mechanical Fermi's Golden Rule, which gives the transition rate from one quantum mechanical eigenstate into a continuum of 
states. 

Consider the coupled crust-core dynamics of section 4. The forced motion of the core Alfven modes due to the acceleration of the 
crust, is 



(Al) 



where a„(?/;) is the displacement of the n-th core Alfven harmonic on the flux-surface ip with frequency ct„, ^(il'jra) is the 
acceleration of the crust at the location where the flux surface ip intersects the crust, and Cn{ip) = {w{ip,x)iS,n) is a coupling 
constant (see Eq. (|49]l). Suppose that we keep the system initially fixed in a position where the crust is displaced with amplitude 
bm,o according to the m-th eigenmode, i.e. ^ = &m,oCm> and the continuum oscillators are relaxed; «„(■(/;) = 0. At time t = 
we release the crust which starts oscillating at frequency ftm- Suppose that the damping timescale „i of the crustal mode is 
much larger than its period = 27r/51m, then the crust oscillates at roughly constant amplitude, i.e. bm{t) ~ 6m, o cos ftmt- This 
motion forces the Alfven oscillators according to 



dn{lp) + Cr^(V')a„(^) = Cn{lp) f^^^m^O Cm ("0, f'o) COS 



(A2) 



One can solve the time-evolution of the oscillator a„(<) using standard techniques (see e.g. Landau & Lifshitz, Mechanics §22). 
After a time t the energy per flux surface Enii') = l/2(a^ + c^n'^n) absorbed by the oscillator is 



COS n,r,t'e-"^"*'dt' 



(A3) 



It is easy to verify that at late times the term between the vertical brackets in Eq. (A3 1 becomes narrowly peaked around (7„ = 



so that the bulk of energy is transported to oscillators which are in (near) resonance with the crust. The average rate of energy (per 
flux surface) transfer (£n('0, t)) from the crust to the to the flux surface at time t is £n{ip, t)/t. For sufficiently large t one finds 



(A4) 



where S{rtm — cr„) is a Dirac delta function. This expression is exact in the limit of t — )• oo. The total rate of energy transfer E 



from the crust to the Alfven continuum is then obtained simply by integrating Eq. ( A4 1 over il> and summing over all n 



(A5) 



here ipk denotes flux surfaces that are in resonance with the crustal motion, <Jn{tpk) = ^m- Since for a given n, the crustal mode 
may be in resonance with Alfven modes in several flux surfaces i/jk, the total energy transfer is obtained by summing over the 



index k. Eq. ( A5 i, which is the analog of the quantum physics' Fermi's Golden Rule, leads to a simple expression for the energy 
damping timescale „j (=1/2 Trf_„, ) of the crustal mode 



Tea 



E{t = 0) 
E 



E 



■4'=ipk 



(A6) 



where E{t = 0) = l/2i7^,j6^„ q is the initial energy of the m-th crustal mode. Using numerical simulations, we verified the 
correctness of Eq. ( A6 1. Even for very short damping times, i.e. Td = 2te ~ 27r/wo, Eq. ( A6 1 proves remarkably accurate. 
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